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Abstract 

An efficient new method is presented to calculate the quantum transports using periodic bound- 
ary conditions. This method allows the use of conventional ground state ab initio programs without 
big changes. The computational effort is only a few times of a normal ground state calculation, 
thus it makes accurate quantum transport calculations for large systems possible. 

PACS numbers: 71.15.-m, 73.63.-b, 73.22.-f 



Quantum transport for molecules, nanowires, and nanodevices is a fast growing research 
area in both experiment and theory, with the potential of replacing the current Si based 
technology after the Moor's law reaches its limit in about 15 years. In the theoretical ballistic 
transport calculations, a key step is to calculate the current via the Landauer formula: 

/=-/ j:Tn{E)dE, (1) 

where /i^ and are left and right electrode Fermi energies (assuming the current flows 
from right to left in z direction), and Tn{E) is the transmission coefficient for the nth right 
hand electrode channel (band) at energy E. There are two major ways to calculate Tn{E). 
One is to use the Green's function G{r, r', E) of the system. However since G is a double 
variable function, computationally this approach can be quite expensive, thus it is mostly 
used for localized basis set methods The other way to calculate Tn{E) is to solve the 
following scattering states: 

HtPUr) = EtPUr) (2) 

and for z oo(— oo): 

i^scir) = E[^f "Vf "^*(r) + B^^^U^^'\r)], (3) 

n 

with conditions: = except = 1 for one m, and = assummg 

dE^^^\k)/dk > 0). In above, H is the single particle Hamiltonian, and (j)n^^\r) = 
Un,k„ii^)^xp{ik^^^^ z) are right going running waves in the the right(R) and left(L) elec- 
trodes, and (pn^^^* the left going running waves. E^^^^^kj^^^^) = E are the elec- 
trode band structure. The transmission coefficient for channel m can be calculated as 



Tm{E) = [J2n\An?{dE^{k)/dk)]/{dE^{k)/dk), and the reflection coefficient can 
lated as Rm{E) = [En \Bn?{dE^{k)Jdk)]/ {dE^{k) / dk). Transfer matrix method 



the Lippmann-Schwinger equation 



3e calcu- 
and 



have been used to solve Eqs(2),(3). Unfortunately, 



the transfer matrix method is rather complicated and computationally expensive to deal 
with the nonlocal pseudopotentials Q] and it is often plagued by the numerical instability 



due to the evanescent states in a multi-channel electrode 



. On the other hand, the use of 



Lippmann-Schwinger equation requires the solution of a linear equation of the dimension 
of the full system, and it also needs the Green's function of the two electrode system un- 



der a potential bias. As a result, currently this approach is only used for jellium electrode 
model and relatively small systems. Overall, compared to the more matured ground state 
calculations, all the current methods for transport calculations are complicated and compu- 
tationally expensive, and they can only be used to calculate relatively small systems although 
there is a strong need to study the transports of large molecules and nanostructures. Here, 
we present a new and simple approach which makes the transport calculation similar to the 
ground state calculation. In this approach, conventional periodic supercell methods and a 
specially designed perturbative approach are used to solve Eqs(2),(3). This allows us to use 
modern ab initio total energy programs without much change. The computational effort is 
similar to a normal ground state local density approximation (LDA) calculation, hence it 
opens the door for quantum transport studies for large systems. 

To demonstrate our method, we have chosen a benzene molecule connected by two Cu 
quantum wires as shown in Fig.l. This system is chosen since the conductivity of the benzene 
molecule is well studied M. M and similar quantum wires have been used as electrodes in 

Q 

previous quantum transport calculations [3]. Two hydrogen atoms at the two ends of the 
benzene molecule are replaced by two sulfur atoms, which are bonded to two central Cu 
atoms at the electrode. The atomic positions of the molecule are relaxed at the zero voltage 
bias under LDA calculation. We have used norm conservingpseudopotentials and 30 Ryd 
planewave cutoff with a standard planewave LDA program j8|. We have included 5 and 6 
unit cells in the left and right electrodes respectively, and they are connected at boundary B 
in Fig. 1 by periodic boundary condition. The x, y dimensions of the supercell are 3 times the 
width of the Cu wire to avoid possible neighbore-neighbore interactions. After the Kohn- 
Sham single particle potential Vo(r) is obtained from a LDA selfconsistent calculation at the 
zero bias, we have added a potential V/2sin{7Tz/ L') in the central region of the molecule and 
shifted the rest of the right (left) electrode by V/2 {—V/2) to get the potential Vv{r) for a 
bias V system. Although a selfconsistent treatment can be achieved straight forwardly under 
the current approach [since the scattering states of Eq(2) will be calculated], the current 
nonselfconsistent treatment for finite bias V is sufficient in illustrating the new methodology. 
Note that, there is a jump of Vv{r) at the boundary B, but that is not a problem in our 
numerical calculations. 

Figure 2 shows the band structures En{kz) of the quantum wire electrode. There are 9 Cu 
atoms in each unit cell of the electrode. To simplify our calculation and analysis, we did not 



include the 3d electrons in our pseudopotential. Although this will introduce a significant 
error in the electrode total energy, the electronic structure near the Fermi surface and the 
related transport properties are intact. However, for bias larger than 2 V, our electrode 
should be considered as a model electrode due to the lack of Cu 3d electrons. 

To solve the scattering states of Eqs(2),(3), we first calculate the eigenstates {ipiir), Ei} 
of the periodic supercell under Vv{r) with a point (e.g, = Tr/2Lz, where is the 
length of the supercell and this is not the of the electrode as in Fig. 2) using our 
standard LDA program Let's first assume that, using some methods, we can generate 
/ degenerated states ipi^{i){r) they all satisfy the Schrodinger's equation Hipi^ii) = Eitpi^ii) 
(however, for our purpose, this equation needs only to be satisfied within the interior of the 
supercell, not near the boundary B of Fig.l). The idea is to use a linear combination of 
these states to construct the scattering states of Eq(2). First, within the R(L) electrode, 
can be decomposed into the electrode states 0^*-^^(r) just as in Eq(3). From a given 
Ei, the available n and fc^^^-* can be found from Fig.2, or say: En{k^^^^) + fJ'R{L) = E^. 
The corresponding (j)^^^^ is then generated by numerical interpolations from pre-calculated 
electrode states. The expansion coefficients A^^^\i, I), B^^^\i, I) for wavefunction can 
be easily calculated from the functional products like: il)i^(i)(j)*^d^r and (pm'Pl^d^^ ■, where 
f2 is one electrode unit cell at the middle of the electrodes as shown in Fig.l. We found that, 
after ^ the po.Me evanescent state, fl. this expansion typically captures .ore than 
99.99% of the weight of the original '?/'«,(/). The next step is to make a linear combination of 
to get the scattering state ipi^sc of Eq(2): 



^i,sc = E/ Cii)^i) (4) 

here the second equation and the R and L are for r within the right and left electrodes 
respectively. Note that, due to the use of supercell point and both ipi and ip* are used 
as ipi^{i), ipi^sc is no longer periodic at the supercell boundary B. To make ipi,sc in Eq(4) 
as one scattering state of Eq(3), we need to make it satisfy the boundary conditions of 
Eq(3):74^_^^ = 0, Bf^ = 0, by selecting C/. In order to have a solution for C/, we need N 
independent ipi,ii) = 1; A^); if there are N nonzero contributing electrode states n in Eq(4) 
[counting both the left and right electrodes, but 0„ and 0* are counted as one]. 



Notice that, since both il>i{r) and il>*{r) satisfy the Schrodinger's equation Hip — Eiip, 
for systems with only a single channel (N=2), the eigen state ipi{r) itself is enough to 
construct the scattering states V'i.sc from Eq(4). So, the main task here is for multi-channel 
cases. In our example, from Fig.2 we see that for a given energy Ei, we could have 4-5 
channels. To get the N degenerated ipi^(i), we will use a perturbative approach. First, we 
will construct Wm{r) — Um,r{r) (the k — mth electrode state) when z is within the last 
unit cell of the right electrode near boundary B of Fig.l, and Wm{r) — for all the other 
z (see Fig.l). We will add P\Wjn >< Wjn\ as a perturbation in the original H, and solve 
the following eigenstate equation using our standard LDA program (e.g, using conjugate 
gradient method): 

HiPl^ + (3< WM^^ >Wr„^ {E, + AE,,mH^^. (5) 

Here /3 is a very small number, hence A£^j ^ and A'0j,m = '^im~'^i both small. Suppose 

we have solved the above equations for two different m's: mi and 1712- Then we can construct 
ipi,(i) = FiAipi^rni + E2AilJi^rn2, with EiAEi^^i + F2AEi^ra2 = 0. After dropping the second 
order terms AEi^rni,2^'4^i,mi,2 have: 

Notice that the Wm^ terms are nonzero only near the boundary B, so for all the other 
places, we have Hip^i^ — Eiipi^(^iy Thus, ipi,{i) are the wavefunctions we needed. In the 
simple cases, when there are N total electrode states 0^^^^ with nonzero components in 
the expansion of ipi{r) in Eq(4), there will be N/2 right electrode states cj)^. Then the 
perturbations by the related N/2 Wm states (which have the same characteristics and cross 
section symmetries as 0^) will introduce N/2 indepedent perturbative wavefunction changes 
Aipi^^. These Aip^^.^ will generate (A^/2 — 1) independent states (besides the original 
ipi). Thus, the total number of ipi,{i) states (counting also ip*Q^) is just N, the exact number 
we need to construct the scattering state ipi^sc from Eq(4) . This argument remains true when 
there are evanescent states or the number of electrode states in the left and right electrodes 
are not the same. Thus, using this procedure, we are guaranteed that there will be enough 
ipi,{i) states for a given ipi to generate a few corresponding scattering states V'j.sc- 

From the supercell eigenstates {ipi{r), E^}, we can generate a set of {k^} from En{k^) + 



fiR = Ei. These {k^} are shown in Fig. 2 as the crosses for a IV bias case (using all the 
ipi{r) with Ei between the two horizontal arrows in Fig. 2). As can be described by a phase 



accumulation model 



id ], on each band, these have roughly equal distances and their total 
number roughly equals the number of electrode unit cells. We have typically used 6 F point 
electrode states as Wm in Eq(5) starting from the lowest band as annotated in Fig.2. This 
means we have to solve {ipi of Eq(5) 6 times using {ipi} as the initial wavefunctions (Notice 
that, this number 6 is roughly the number of channels in the problem. The same prefactor is 
needed in the calculations of other methods like the transfer matrix or Lippmann-Schwinger 
equation). After {ip^rn} calculated, using Eq(4), we can construct a scattering state 
from each of these k^ shown in Fig.2. Two of these constructed scattering states are shown 
in Fig. 3. Notice that the dash lines are Z); C;'?/'j,(0 of Eq(4), while the solid lines are the 
electrode state decompositions [the second line of the Eq(4)]. Within the electrode, the 
electrode state decomposition gives a very accurate description of the total wavefunction. 
From these scattering states, the transmission coefficients Tm{Ei) can be calculated, and are 
shown in Fig.4 as the symbols. The calculated Tm{Ei) + Rm{Ei) is typically very close to 1, 
indicating the numerical stability of the current method. 

Notice that, unlike the other approaches discussed before where the scattering states 
of arbitrary energy E can be solved, here only the scattering states of energy {Ei} are 
calculated. This translates into finite number of {k^} points as shown in Fig.2 and Fig.4. In 
Eq(l), we need all the energies between fi^ and fiR. This is an complete analogy with the k- 
point integration problem in conventional ground state bulk calculation [the energy integral 
in Eq(l) can also be changed into a k-point integral of the electrode band structure of Fig.2]. 
Thus, similar to the conventional bulk calculations, here we will use an interpolation scheme 
to carry out the integral in Eq(l). First, if the number of {k^} is not enough in Fig.4, we 
can choose a different supercell in our supercell calculations (or change the potential 
Vv{r) near boundary B), and repeat the above procedure. That will give us more {Ei} and 
{kn} points. In our system, we find one calculation is sufficient. We have used a smooth 
curve fn{k) to interpolate the ln{Tn{k^)) points shown in Fig.4. More specifically, we have 
minimized: \ln{Tn{k^)) - fn{k^)\^ + ^ ^ \<P fn{k) / dk'^\'^dk, with 7 ~ 1. Numerically, 
this corresponds to a simple linear equation with discretized k points. The resulting curve is 
shown in Fig.4 for the IV bias case. Using these fn{k), we can calculate the total transmission 
T{E) = J2n'^n{E) of the system. The results are shown in Fig. 5 for different biases. We see 



that T{E) is influenced strongly by two factors. One is the relative energy levels between the 
electrode states and the molecule states. When the bias increases, the molecular levels drop 
relative to the right electrode state levels. As a result, the magnitude of the transmission 
decreases near the region of -3 eV. Another factor is the band structure of the electrode. 
There is a well shape of T{E) near -0.3 eV. This is caused by band gaps of the 2,3 bands at 
the X' point in Fig. 2. The T{E) also shows a big drop at -3.3 eV. This is due to the end of 
2,3 bands at the F point. 

After T{E)'s of Fig. 5 are obtained, a simple energy integration between -V to will give 
us the total current /. The resulting I{V), and the conductance dl/dV are shown in Fig.6. 
We do see the well known peak and dip for this system in the conductance around 2V. Our 
calculated peak and dip positions of 1.8V and 2.3V corresponds well with the experimental 
results of 1.4V and 2.4V |6:], and this agreement is better than previously calculated results 
P]. We also see the marks of the electrode electronic structures. Near 0.3V, the conductance 
shows a well shape, again due to the band gaps of 2,3 bands at the X' point. Above 3.3V, we 
see a big drop, then the negative conductance. The drop is due to the end of the 2,3 bands 
at the F points, and the negative conductance is because the conducting electrode levels 
(energy window) are moving away from the conducting molecular levels. Here we see that, 
the electronic structure of the electrode is extremely important in determining the overall 
conductance of the system. 

In summary, we have presented a simple and numerically stable scheme to calculate the 
quantum transport. Within this scheme, the conventional periodic supercells and ground 
state ab initio programs can be used without much change. The computational effort is only 
a few times of a conventional ground state calculation. This promised quantum transport 
calculations for much larger systems which cannot be tackled by othe methods. The imple- 
mentation of this method is simple and straight forward based on any conventional ground 
state ab initio programs. 
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FIG. 1: A schematic view of the calculated system. 



FIG. 2: The band structure of the electrode. Each continuous line from F to X' is denoted as one 
band. The zero is the electrode Fermi energy. The crosses are the points, see text for details. 

FIG. 3: The constructed scattering states from Eq(4). The dashed and solid lines correspond to 
the first and second lines of Eq(4) respectively. T is the transmission coefficient, E is the eigen 
energy, and band number indicate the n of (/in in Eq(4). The bias of the system is IV. 



FIG. 4: The calculated transmission coefficients T„(A;^) (symbols) and fitted smooth curves fn{k) 
(lines) for different bands of the electrode. The bias of the system is IV. 



FIG. 5: The calculated total transmission coefficients T{E) (which can be larger than 1) for 
different biases. The zero is the right electrode Fermi energy. For a given bias V, there are net 
right to left current flow only within the [—V, 0] energy window. 



FIG. 6: The calculated I-V curve and the corresponding conductance. Go = 2e/h = 77/xs. 
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